### load data

library (tidyverse)
library (ggplot2)
library(jtools)
library (officer)
library (flextable)
library (broom.mixed)
library (car)
library(dplyr)
library (ggplot2)
library (ggeffects)
library(cowplot)
library (psych)
library(readr)
library(lmtest)
library(sandwich)

df_clean <- na.omit(data)

df_clean$social_media_binary <- as.factor(df_clean$social_media_binary)
df_clean$express_binary <- as.factor(df_clean$express_binary)
df_clean$gender <- as.factor(df_clean$gender)
df_clean$employed <- as.factor(df_clean$employed)
df_clean$education <- as.factor(df_clean$education)
df_clean$urban <- as.factor(df_clean$urban)
df_clean$country_name <- as.factor(df_clean$country_name)
df_clean$wave <- as.factor(df_clean$wave)


# Political Trust Model
model_political <- lm(political_trust_index ~ social_media_binary + express_binary + age + gender +
                        employed + education + income + urban + country_name + wave,
                      data = df_clean)

# Protective Trust Model
model_protective <- lm(protective_trust_index ~ social_media_binary + express_binary + age + gender +
                         employed + education + income + urban + country_name + wave,
                       data = df_clean)


# Robust standard errors clustered by country
political_clustered <- coeftest(model_political, vcov = vcovCL(model_political, cluster = ~country_name))
protective_clustered <- coeftest(model_protective, vcov = vcovCL(model_protective, cluster = ~country_name))


print(political_clustered)

print(protective_clustered)

# Model Diagnostics

vif(model_political)
vif(model_protective)
# no multicolinnearity 

# political trust model
plot(model_political$fitted.values, residuals(model_political),
     main = "Political Trust Model: Residuals vs Fitted",
     xlab = "Fitted values", ylab = "Residuals")
abline(h = 0, col = "red")

# protective trust model
plot(model_protective$fitted.values, residuals(model_protective),
     main = "Protective Trust Model: Residuals vs Fitted",
     xlab = "Fitted values", ylab = "Residuals")
abline(h = 0, col = "red")


qqnorm(residuals(model_political), main = "Political Trust Model: Q-Q Plot")
qqline(residuals(model_political), col = "red")

qqnorm(residuals(model_protective), main = "Protective Trust Model: Q-Q Plot")
qqline(residuals(model_protective), col = "red")


bptest(model_political)   
bptest(model_protective)

# tables and plots

library(modelsummary)


modelsummary(
  list("Political Trust" = model_political, "Protective Trust" = model_protective),
  vcov = list(vcovCL(model_political, cluster = ~country_name),
              vcovCL(model_protective, cluster = ~country_name)),
  statistic = "conf.int",
  stars = TRUE,
  gof_omit = "AIC|BIC|Log.Lik",
  output = "trust_models.docx"
)



# ORDINAL MODEL

library(ordinal)


df_clean$ngos <- factor(df_clean$ngos, levels = c(1, 2, 3, 4), ordered = TRUE)


model_ngos <- clm(ngos ~ social_media_binary + express_binary + age + gender +
                    employed + education + income + urban + country_name + wave,
                  data = df_clean)


library(clubSandwich)


vcov_ngos <- vcovCL(model_ngos, cluster = df_clean$country_name)
coeftest(model_ngos, vcov = vcov_ngos)


modelsummary(
  list("Trust in NGOs" = model_ngos),
  vcov = list(vcovCL(model_ngos, cluster = ~country_name)),
  statistic = "conf.int",
  stars = TRUE,
  gof_omit = "AIC|BIC|Log.Lik",
  output = "ngo_trust_model.docx"
)


modelsummary(
  model_ngos,
  statistic = "conf.int",
  stars = TRUE,
  gof_omit = "AIC|BIC|Log.Lik",
  coef_omit = "Intercept|\\|",  
  output = "trust_ngos_clean.docx"
)


# plots

library(ggplot2)
library(dplyr)


political_df <- data.frame(
  variable = c("social_media_binary1", "express_binary1", "age", "gender1", "employed1",
               "education2", "education3", "income", "urban1", "country_nameIndonesia",
               "country_nameMalaysia", "country_nameMyanmar", "country_namePhilippines",
               "country_nameSingapore", "country_nameThailand", "country_nameVietnam", "wave5"),
  estimate = c(-0.005, 0.470, -0.0049, -0.167, 0.034,
               0.865, 1.538, 0.0005, -0.151, -0.931,
               -0.331, -0.204, 0.664, -1.665, 0.238,
               -2.825, -2.497),
  std_error = c(0.198, 0.122, 0.0047, 0.0942, 0.133,
                0.455, 0.486, 0.0831, 0.252, 0.380,
                0.381, 0.359, 0.333, 0.325, 0.338,
                0.388, 0.667)
)


political_df <- political_df %>%
  mutate(
    lower = estimate - 1.96 * std_error,
    upper = estimate + 1.96 * std_error,
    label = factor(c("Social Media Use", "Political Expression", "Age", "Gender (Female)", "Employment Status",
                     "Education (Secondary)", "Education (College)", "Income", "Urban", "Indonesia",
                     "Malaysia", "Myanmar", "Philippines", "Singapore", "Thailand", "Vietnam", "Wave 5"),
                   levels = rev(c("Social Media Use", "Political Expression", "Age", "Gender (Female)", "Employment Status",
                                  "Education (Secondary)", "Education (College)", "Income", "Urban", "Indonesia",
                                  "Malaysia", "Myanmar", "Philippines", "Singapore", "Thailand", "Vietnam", "Wave 5")))
  )

# Plot
ggplot(political_df, aes(x = estimate, y = label)) +
  geom_point(color = "black") +
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2, color = "darkgray") +
  geom_text(aes(label = round(estimate, 2)), 
            vjust = -1.2, hjust = 0.5, size = 3, color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "Political Trust Model Coefficients",
       x = "Coefficient Estimate", y = NULL) +
  theme_minimal()



library(ggplot2)
library(dplyr)

# Data for protective trust model
protective_df <- data.frame(
  variable = c("social_media_binary1", "express_binary1", "age", "gender1", "employed1",
               "education2", "education3", "income", "urban1", "country_nameIndonesia",
               "country_nameMalaysia", "country_nameMyanmar", "country_namePhilippines",
               "country_nameSingapore", "country_nameThailand", "country_nameVietnam", "wave5"),
  estimate = c(0.015, 0.084, -0.0016, 0.011, 0.074,
               0.206, 0.374, 0.040, -0.083, -0.640,
               -0.412, 0.611, 0.137, -0.626, -0.074,
               -0.872, -0.951),
  std_error = c(0.063, 0.068, 0.002, 0.042, 0.065,
                0.077, 0.088, 0.031, 0.092, 0.138,
                0.123, 0.119, 0.099, 0.096, 0.110,
                0.127, 0.200)
)


protective_df <- protective_df %>%
  mutate(
    lower = estimate - 1.96 * std_error,
    upper = estimate + 1.96 * std_error,
    label = factor(c("Social Media Use", "Political Expression", "Age", "Gender (Female)", "Employment Status",
                     "Education (Secondary)", "Education (College)", "Income", "Urban", "Indonesia",
                     "Malaysia", "Myanmar", "Philippines", "Singapore", "Thailand", "Vietnam", "Wave 5"),
                   levels = rev(c("Social Media Use", "Political Expression", "Age", "Gender (Female)", "Employment Status",
                                  "Education (Secondary)", "Education (College)", "Income", "Urban", "Indonesia",
                                  "Malaysia", "Myanmar", "Philippines", "Singapore", "Thailand", "Vietnam", "Wave 5")))
  )

# Plot
ggplot(protective_df, aes(x = estimate, y = label)) +
  geom_point(color = "black") +
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2, color = "darkgray") +
  geom_text(aes(label = round(estimate, 2)), 
            vjust = -1.2, hjust = 0.5, size = 3, color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "Protective Trust Model Coefficients",
       x = "Coefficient Estimate", y = NULL) +
  theme_minimal()


library(ggplot2)
library(dplyr)
library(patchwork)  

labels <- c("Social Media Use", "Political Expression", "Age", "Gender (Female)", "Employment Status",
            "Education (Secondary)", "Education (College)", "Income", "Urban", "Indonesia",
            "Malaysia", "Myanmar", "Philippines", "Singapore", "Thailand", "Vietnam", "Wave 5")

# Political model
political_df <- data.frame(
  label = factor(labels, levels = rev(labels)),
  estimate = c(-0.005, 0.470, -0.0049, -0.167, 0.034,
               0.865, 1.538, 0.0005, -0.151, -0.931,
               -0.331, -0.204, 0.664, -1.665, 0.238,
               -2.825, -2.497),
  std_error = c(0.198, 0.122, 0.0047, 0.0942, 0.133,
                0.455, 0.486, 0.0831, 0.252, 0.380,
                0.381, 0.359, 0.333, 0.325, 0.338,
                0.388, 0.667)
) %>%
  mutate(lower = estimate - 1.96 * std_error,
         upper = estimate + 1.96 * std_error)

# Protective model
protective_df <- data.frame(
  label = factor(labels, levels = rev(labels)),
  estimate = c(0.015, 0.084, -0.0016, 0.011, 0.074,
               0.206, 0.374, 0.040, -0.083, -0.640,
               -0.412, 0.611, 0.137, -0.626, -0.074,
               -0.872, -0.951),
  std_error = c(0.063, 0.068, 0.002, 0.042, 0.065,
                0.077, 0.088, 0.031, 0.092, 0.138,
                0.123, 0.119, 0.099, 0.096, 0.110,
                0.127, 0.200)
) %>%
  mutate(lower = estimate - 1.96 * std_error,
         upper = estimate + 1.96 * std_error)


plot_model <- function(df, title) {
  ggplot(df, aes(x = estimate, y = label)) +
    geom_point(color = "black") +
    geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2, color = "darkgray") +
    geom_text(aes(label = round(estimate, 2)), vjust = -1.2, size = 3, color = "black") +
    geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
    labs(title = title, x = "Coefficient Estimate", y = NULL) +
    theme_minimal()
}


plot1 <- plot_model(political_df, "Political Trust Model")
plot2 <- plot_model(protective_df, "Protective Trust Model")


combined_plot <- plot1 + plot2
combined_plot


# plot ordinal model


library(ggplot2)
library(dplyr)

# Original data
ordinal_df <- data.frame(
  variable = c("social_media_binary1", "express_binary1", "age", "gender1", "employed1",
               "education2", "education3", "income", "urban1", "country_nameIndonesia",
               "country_nameMalaysia", "country_nameMyanmar", "country_namePhilippines",
               "country_nameSingapore", "country_nameThailand", "country_nameVietnam", "wave5"),
  estimate = c(-0.036, -0.127, -0.0012, 0.027, 0.056,
               0.154, 0.298, 0.018, -0.051, 0.125,
               0.913, 0.638, 1.024, 1.015, 1.261,
               0.482, -1.235),
  std_error = c(0.0593, 0.0494, 0.0023, 0.0481, 0.0505,
                0.1089, 0.1344, 0.0221, 0.1026, 0.1320,
                0.1373, 0.1190, 0.1195, 0.1223, 0.1143,
                0.1263, 0.1977)
)


label_map <- c(
  social_media_binary1 = "Social Media Use",
  express_binary1 = "Political Expression",
  age = "Age",
  gender1 = "Gender (Female)",
  employed1 = "Employment Status",
  education2 = "Education (Secondary)",
  education3 = "Education (College)",
  income = "Income",
  urban1 = "Urban",
  country_nameIndonesia = "Indonesia",
  country_nameMalaysia = "Malaysia",
  country_nameMyanmar = "Myanmar",
  country_namePhilippines = "Philippines",
  country_nameSingapore = "Singapore",
  country_nameThailand = "Thailand",
  country_nameVietnam = "Vietnam",
  wave5 = "Wave 5"
)


ordinal_df <- ordinal_df %>%
  mutate(
    label = factor(label_map[variable], levels = rev(label_map[variable])),
    lower = estimate - 1.96 * std_error,
    upper = estimate + 1.96 * std_error
  )

# Plot
ggplot(ordinal_df, aes(x = estimate, y = label)) +
  geom_point(color = "black") +
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2, color = "darkgray") +
  geom_text(aes(label = round(estimate, 2)), 
            vjust = -1.2, size = 3, color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "Trust in NGOs",
       x = "Estimate", y = NULL) +
  theme_minimal()



################# REVERSED CODED #############################

# Political Trust Model
model_political_rev <- lm(political_trust_index_rev ~ social_media_binary + express_binary + age + gender +
                        employed + education + income + urban + country_name + wave,
                      data = df_clean)

# Protective Trust Model
model_protective_rev <- lm(protective_trust_index_rev ~ social_media_binary + express_binary + age + gender +
                         employed + education + income + urban + country_name + wave,
                       data = df_clean)


# Robust standard errors clustered by country
political_clustered <- coeftest(model_political_rev, vcov = vcovCL(model_political_rev, cluster = ~country_name))
protective_clustered <- coeftest(model_protective_rev, vcov = vcovCL(model_protective_rev, cluster = ~country_name))


print(political_clustered)

print(protective_clustered)

# Model Diagnostics

vif(model_political_rev)
vif(model_protective_rev)
# no multicolinnearity 

library(modelsummary)


modelsummary(
  list("Political Trust" = model_political_rev, "Protective Trust" = model_protective_rev),
  vcov = list(vcovCL(model_political_rev, cluster = ~country_name),
              vcovCL(model_protective_rev, cluster = ~country_name)),
  statistic = "conf.int",
  stars = TRUE,
  gof_omit = "AIC|BIC|Log.Lik",
  output = "trust_models_REV.docx"
)


# ORDINAL MODEL REVERSED

library(ordinal)

# Convert ngo_trust to ordered factor (1 = trust, 4 = no trust)
df_clean$ngos_rev <- factor(df_clean$ngos_rev, levels = c(1, 2, 3, 4), ordered = TRUE)

# Fit the ordinal logistic regression model
model_ngos_rev <- clm(ngos_rev ~ social_media_binary + express_binary + age + gender +
                    employed + education + income + urban + country_name + wave,
                  data = df_clean)


library(clubSandwich)

# Clustered robust standard errors
vcov_ngos <- vcovCL(model_ngos_rev, cluster = df_clean$country_name)
coeftest(model_ngos_rev, vcov = vcov_ngos)


modelsummary(
  list("Trust in NGOs" = model_ngos_rev),
  vcov = list(vcovCL(model_ngos_rev, cluster = ~country_name)),
  statistic = "conf.int",
  stars = TRUE,
  gof_omit = "AIC|BIC|Log.Lik",
  coef_omit = "Intercept|\\|",
  output = "ngo_trust_model_rev.docx"
)


# PLOTTING

library(ggplot2)
library(dplyr)
library(patchwork)  # for combining plots

# Data for Political Trust Model
political_df <- data.frame(
  variable = c("social_media_binary1", "express_binary1", "age", "gender1", "employed1",
               "education2", "education3", "income", "urban1", "country_nameIndonesia",
               "country_nameMalaysia", "country_nameMyanmar", "country_namePhilippines",
               "country_nameSingapore", "country_nameThailand", "country_nameVietnam", "wave5"),
  estimate = c(0.0053, -0.4705, 0.0049, 0.1667, -0.0339,
               -0.8647, -1.5385, -0.0005, 0.1510, 0.9307,
               0.3314, 0.2039, -0.6635, 1.6652, -0.2383,
               2.8250, 2.4974),
  std_error = c(0.1977, 0.1220, 0.0047, 0.0942, 0.1331,
                0.4554, 0.4862, 0.0831, 0.2519, 0.3803,
                0.3811, 0.3593, 0.3326, 0.3252, 0.3378,
                0.3882, 0.6670)
)

# Data for Protective Trust Model
protective_df <- data.frame(
  variable = c("social_media_binary1", "express_binary1", "age", "gender1", "employed1",
               "education2", "education3", "income", "urban1", "country_nameIndonesia",
               "country_nameMalaysia", "country_nameMyanmar", "country_namePhilippines",
               "country_nameSingapore", "country_nameThailand", "country_nameVietnam", "wave5"),
  estimate = c(-0.0147, -0.0840, 0.0016, -0.0112, -0.0745,
               -0.2059, -0.3740, -0.0399, 0.0826, 0.6405,
               0.4116, -0.6106, -0.1375, 0.6256, 0.0737,
               0.8722, 0.9505),
  std_error = c(0.0631, 0.0683, 0.0020, 0.0421, 0.0645,
                0.0766, 0.0878, 0.0312, 0.0925, 0.1377,
                0.1230, 0.1195, 0.0992, 0.0959, 0.1100,
                0.1266, 0.1996)
)

# Function to prepare plot
plot_model <- function(df, title) {
  df <- df %>%
    mutate(
      lower = estimate - 1.96 * std_error,
      upper = estimate + 1.96 * std_error,
      label = factor(variable, levels = rev(variable))
    )
  
  ggplot(df, aes(x = estimate, y = label)) +
    geom_point(color = "blue") +
    geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2, color = "gray") +
    geom_text(aes(label = round(estimate, 2)), vjust = -1.2, size = 3, color = "black") +
    geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
    labs(title = title, x = "Coefficient Estimate", y = NULL) +
    theme_minimal()
}

# Create plots
plot1 <- plot_model(political_df, "Political Trust Model (Reversed Scale)")
plot2 <- plot_model(protective_df, "Protective Trust Model (Reversed Scale)")

# Combine plots side by side
combined_plot <- plot1 + plot2
combined_plot


library(ggplot2)
library(dplyr)
library(patchwork)

# Define descriptive labels
label_map <- c(
  social_media_binary1 = "Social Media Use",
  express_binary1 = "Political Expression",
  age = "Age",
  gender1 = "Gender (Female)",
  employed1 = "Employment Status",
  education2 = "Education (Secondary)",
  education3 = "Education (College)",
  income = "Income",
  urban1 = "Urban",
  country_nameIndonesia = "Indonesia",
  country_nameMalaysia = "Malaysia",
  country_nameMyanmar = "Myanmar",
  country_namePhilippines = "Philippines",
  country_nameSingapore = "Singapore",
  country_nameThailand = "Thailand",
  country_nameVietnam = "Vietnam",
  wave5 = "Wave 5"
)

# Apply labels to both data frames
political_df <- political_df %>%
  mutate(
    label = factor(label_map[variable], levels = rev(label_map[variable])),
    lower = estimate - 1.96 * std_error,
    upper = estimate + 1.96 * std_error
  )

protective_df <- protective_df %>%
  mutate(
    label = factor(label_map[variable], levels = rev(label_map[variable])),
    lower = estimate - 1.96 * std_error,
    upper = estimate + 1.96 * std_error
  )

# Plot function
plot_model <- function(df, title) {
  ggplot(df, aes(x = estimate, y = label)) +
    geom_point(color = "black") +
    geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2, color = "darkgray") +
    geom_text(aes(label = round(estimate, 2)), vjust = -1.2, size = 3, color = "black") +
    geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
    labs(title = title, x = "Coefficient Estimate", y = NULL) +
    theme_minimal()
}

# Create and combine plots
plot1 <- plot_model(political_df, "Political Trust Model")
plot2 <- plot_model(protective_df, "Protective Trust Model")
combined_plot <- plot1 + plot2
combined_plot


library(ggplot2)
library(dplyr)


ngo_df <- data.frame(
  variable = c("social_media_binary1", "express_binary1", "age", "gender1", "employed1",
               "education2", "education3", "income", "urban1", "country_nameIndonesia",
               "country_nameMalaysia", "country_nameMyanmar", "country_namePhilippines",
               "country_nameSingapore", "country_nameThailand", "country_nameVietnam", "wave5"),
  estimate = c(0.0362, 0.1274, 0.0012, -0.0266, -0.0561,
               -0.1541, -0.2982, -0.0182, 0.0510, -0.1248,
               -0.9133, -0.6377, -1.0238, -1.0151, -1.2615,
               -0.4819, 1.2345),
  std_error = c(0.0593, 0.0494, 0.0023, 0.0481, 0.0505,
                0.1089, 0.1344, 0.0221, 0.1026, 0.1320,
                0.1373, 0.1190, 0.1195, 0.1223, 0.1143,
                0.1263, 0.1977)
)


ngo_df <- ngo_df %>%
  mutate(
    lower = estimate - 1.96 * std_error,
    upper = estimate + 1.96 * std_error
  )

label_map <- c(
  social_media_binary1 = "Social Media Use",
  express_binary1 = "Political Expression",
  age = "Age",
  gender1 = "Gender (Female)",
  employed1 = "Employment Status",
  education2 = "Education (Secondary)",
  education3 = "Education (College)",
  income = "Income",
  urban1 = "Urban",
  country_nameIndonesia = "Indonesia",
  country_nameMalaysia = "Malaysia",
  country_nameMyanmar = "Myanmar",
  country_namePhilippines = "Philippines",
  country_nameSingapore = "Singapore",
  country_nameThailand = "Thailand",
  country_nameVietnam = "Vietnam",
  wave5 = "Wave 5"
)


ngo_df <- ngo_df %>%
  mutate(label = factor(label_map[variable], levels = rev(label_map[variable])))

# Plot
ggplot(ngo_df, aes(x = estimate, y = label)) +
  geom_point(color = "black") +
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2, color = "darkgray") +
  geom_text(aes(label = round(estimate, 2)), vjust = -1.2, size = 3, color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "Trust in NGOs",
       x = "Coefficient Estimate", y = NULL) +
  theme_minimal()



############# robustness

library(ordinal)


df_clean$neighbors_rev <- factor(df_clean$neighbors_rev, levels = c(1, 2, 3, 4), ordered = TRUE)


model_neigh_rev <- clm(neighbors_rev ~ social_media_binary + express_binary + age + gender +
                        employed + education + income + urban + country_name + wave,
                      data = df_clean)


library(clubSandwich)


vcov_neigh <- vcovCL(model_neigh_rev, cluster = df_clean$country_name)
coeftest(model_neigh_rev, vcov = vcov_neigh)


modelsummary(
  list("Trust in Neighbors" = model_neigh_rev),
  vcov = list(vcovCL(model_ngos_rev, cluster = ~country_name)),
  statistic = "conf.int",
  stars = TRUE,
  gof_omit = "AIC|BIC|Log.Lik",
  coef_omit = "Intercept|\\|",
  output = "neighbor_trust_model_rev.docx"
)

library(ordinal)


df_clean$others_rev <- factor(df_clean$others_rev, levels = c(1, 2, 3, 4), ordered = TRUE)


model_others_rev <- clm(others_rev ~ social_media_binary + express_binary + age + gender +
                        employed + education + income + urban + country_name + wave,
                      data = df_clean)


library(clubSandwich)


vcov_other <- vcovCL(model_others_rev, cluster = df_clean$country_name)
coeftest(model_others_rev, vcov = vcov_other)


modelsummary(
  list("Trust in Others" = model_others_rev),
  vcov = list(vcovCL(model_ngos_rev, cluster = ~country_name)),
  statistic = "conf.int",
  stars = TRUE,
  gof_omit = "AIC|BIC|Log.Lik",
  coef_omit = "Intercept|\\|",
  output = "others_trust_model_rev.docx"
)

library(ordinal)


df_clean$newspaper_rev <- factor(df_clean$newspaper_rev, levels = c(1, 2, 3, 4), ordered = TRUE)


model_news_rev <- clm(newspaper_rev ~ social_media_binary + express_binary + age + gender +
                        employed + education + income + urban + country_name + wave,
                      data = df_clean)


library(clubSandwich)


vcov_news <- vcovCL(model_news_rev, cluster = df_clean$country_name)
coeftest(model_news_rev, vcov = vcov_news)


modelsummary(
  list("Trust in Newspapers" = model_news_rev),
  vcov = list(vcovCL(model_ngos_rev, cluster = ~country_name)),
  statistic = "conf.int",
  stars = TRUE,
  gof_omit = "AIC|BIC|Log.Lik",
  coef_omit = "Intercept|\\|",
  output = "news_trust_model_rev.docx"
)


library(ordinal)


df_clean$tv_rev <- factor(df_clean$tv_rev, levels = c(1, 2, 3, 4), ordered = TRUE)


model_tv_rev <- clm(tv_rev ~ social_media_binary + express_binary + age + gender +
                        employed + education + income + urban + country_name + wave,
                      data = df_clean)


library(clubSandwich)
library (modelsummary)


vcov_tv <- vcovCL(model_tv_rev, cluster = df_clean$country_name)
coeftest(model_tv_rev, vcov = vcov_tv)


modelsummary(
  list("Trust in Television" = model_tv_rev),
  vcov = list(vcovCL(model_ngos_rev, cluster = ~country_name)),
  statistic = "conf.int",
  stars = TRUE,
  gof_omit = "AIC|BIC|Log.Lik",
  coef_omit = "Intercept|\\|",
  output = "tv_trust_model_rev.docx"
)


library(ordinal)
library(clubSandwich)
library(modelsummary)


vcov_list <- list(
  "Trust in Neighbors" = vcovCL(model_neigh_rev, cluster = ~country_name),
  "Trust in Others" = vcovCL(model_others_rev, cluster = ~country_name),
  "Trust in Newspapers" = vcovCL(model_news_rev, cluster = ~country_name),
  "Trust in Television" = vcovCL(model_tv_rev, cluster = ~country_name)
)


model_list <- list(
  "Trust in Neighbors" = model_neigh_rev,
  "Trust in Others" = model_others_rev,
  "Trust in Newspapers" = model_news_rev,
  "Trust in Television" = model_tv_rev
)

# Create a single Word 
modelsummary(
  models = model_list,
  vcov = vcov_list,
  statistic = "conf.int",
  stars = TRUE,
  coef_omit = "Intercept|\\|",  # remove intercepts and cut points
  gof_omit = "AIC|BIC|Log.Lik",
  output = "trust_robust_models_rev.docx"
)


df_clean$election_rev <- factor(df_clean$election_rev, levels = c(1, 2, 3, 4), ordered = TRUE)

model_elec_rev <- clm(election_rev ~ social_media_binary + express_binary + age + gender +
                      employed + education + income + urban + country_name + wave,
                    data = df_clean)



vcov_elec <- vcovCL(model_elec_rev, cluster = df_clean$country_name)
coeftest(model_elec_rev, vcov = vcov_elec)



df_clean$civil_rev <- factor(df_clean$civil_rev, levels = c(1, 2, 3, 4), ordered = TRUE)


model_civil_rev <- clm(civil_rev ~ social_media_binary + express_binary + age + gender +
                        employed + education + income + urban + country_name + wave,
                      data = df_clean)


vcov_civil <- vcovCL(model_civil_rev, cluster = df_clean$country_name)
coeftest(model_civil_rev, vcov = vcov_civil)



models <- list(
  "Trust in Elections" = model_elec_rev,
  "Trust in Civil Service" = model_civil_rev
)

vcovs <- list(
  "Trust in Elections" = vcov_elec,
  "Trust in Civil Service" = vcov_civil
)


modelsummary(
  models,
  vcov = vcovs,
  statistic = "conf.int",  
  fmt = 3,
  stars = TRUE,
  output = "robust_gov_rev.docx"  
)


















